% 读取Excel数据
data = readtable('附件1 数据-C题.xlsx');

% 数据清洗
% 1. 移除失业注销时间或失业时间为空的记录
data = data(~ismissing(data.c_ajc090), :); % 至少需要失业时间

% 2. 移除关键字段为空的记录
key_fields = {'age', 'sex', 'edu_level', 'c_aac183', 'b_aab022'};
data = data(~any(ismissing(data(:, key_fields)), 2), :);

% 3. 验证年龄、性别、学历的有效性
data = data(data.age >= 0 & data.age <= 100, :); % 年龄在0-100之间
data = data(ismember(data.sex, [1, 2]), :); % 性别为1或2
valid_edu = [10,11,12,13,14,15,16,17,18,19,20,21,22,23,28,30,31,32,33,...
             40,41,42,43,44,45,46,47,48,49,50,60,61,62,63,70,71,73,...
             80,81,83,90,91,92,93,99];
data = data(ismember(data.edu_level, valid_edu), :);

% 4. 转换日期格式并验证
try
    data.c_acc028 = datetime(data.c_acc028, 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
    data.c_ajc090 = datetime(data.c_ajc090, 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
catch
    % 如果日期格式不一致，尝试无时间部分
    data.c_acc028 = datetime(data.c_acc028, 'InputFormat', 'yyyy-MM-dd');
    data.c_ajc090 = datetime(data.c_ajc090, 'InputFormat', 'yyyy-MM-dd');
end
% 移除无效日期
data = data(~isnat(data.c_ajc090), :);

% 确定就业状态
data.employment_status = categorical(~ismissing(data.c_acc028) & (data.c_acc028 > data.c_ajc090), [true, false], {'就业', '失业'});

% 检查是否还有数据
if isempty(data)
    error('清洗后没有有效数据！');
end

% 处理年龄：分组为年龄段
age_bins = [0, 25, 35, 45, Inf];
age_labels = {'<25', '25-35', '35-45', '>45'};
data.age_group = discretize(data.age, age_bins, 'categorical', age_labels);

% 处理性别：转换为描述
sex_labels = {'男', '女'};
data.sex = categorical(data.sex, [1, 2], sex_labels);

% 处理学历：根据编码转换为描述
edu_map = containers.Map(...
    {'10','11','12','13','14','15','16','17','18','19', ... % 研究生及以上
     '20','21','22','23','28', ... % 本科
     '30','31','32','33', ... % 专科
     '40','41','42','43','44','45','46','47','48','49','91','92','93', ... % 中职
     '50','60','61','62','63','70','71','73','80','81','83','90', ... % 高中及以下
     '99'}, ... % 其他
    {'研究生及以上','研究生及以上','研究生及以上','研究生及以上','研究生及以上',...
     '研究生及以上','研究生及以上','研究生及以上','研究生及以上','研究生及以上',...
     '本科','本科','本科','本科','本科',...
     '专科','专科','专科','专科',...
     '中职','中职','中职','中职','中职','中职','中职','中职','中职','中职','中职','中职','中职',...
     '高中及以下','高中及以下','高中及以下','高中及以下','高中及以下',...
     '高中及以下','高中及以下','高中及以下','高中及以下','高中及以下',...
     '高中及以下','高中及以下',...
     '其他'});
data.edu = cell(height(data), 1);
for i = 1:height(data)
    edu_key = num2str(data.edu_level(i));
    if isKey(edu_map, edu_key)
        data.edu{i} = edu_map(edu_key);
    else
        data.edu{i} = '其他';
    end
end
data.edu = categorical(data.edu);

% 处理专业：简化为主要类别
data.prof = data.c_aac183;
data.prof = strrep(data.prof, '其他学科', '其他');
prof_categories = unique(data.prof);
if length(prof_categories) > 5
    prof_counts = groupcounts(data, 'prof');
    [~, idx] = sort([prof_counts.GroupCount], 'descend');
    top_profs = prof_categories(idx(1:min(4, length(idx))));
    data.prof(~ismember(data.prof, top_profs)) = {'其他'};
end
data.prof = categorical(data.prof);

% 处理行业：归并为五个类别
agri_codes = {'A'};
manu_codes = {'C'};
const_mining_codes = {'B', 'E'};
trade_trans_codes = {'F', 'G'};
service_codes = {'D', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'Z'};
data.industry = cell(height(data), 1);
for i = 1:height(data)
    code = string(data.b_aab022{i});
    code_prefix = extractBefore(code, 2);
    if ismember(code_prefix, agri_codes)
        data.industry{i} = '农业及相关产业';
    elseif ismember(code_prefix, manu_codes)
        data.industry{i} = '制造业';
    elseif ismember(code_prefix, const_mining_codes)
        data.industry{i} = '建筑及采矿业';
    elseif ismember(code_prefix, trade_trans_codes)
        data.industry{i} = '批发 批发零售及运输业';
    elseif ismember(code_prefix, service_codes)
        data.industry{i} = '服务业及其他';
    else
        data.industry{i} = '服务业及其他';
    end
end
data.industry = categorical(data.industry);


%% === Step 3. 特征数值化 + 建模准备 ===
data.employment_status_numeric = double(data.employment_status == '就业');
data_table = table(data.age, double(data.sex), double(data.edu), ...
    double(data.prof), double(data.industry), data.employment_status_numeric, ...
    'VariableNames', {'age','sex','edu','prof','industry','employment_status'});

% 划分训练集/测试集
cv = cvpartition(height(data_table), 'HoldOut', 0.3);
train_data = data_table(training(cv), :);
test_data = data_table(test(cv), :);

X_train = train_data{:,1:5};
y_train = train_data.employment_status;
X_test = test_data{:,1:5};
y_test = test_data.employment_status;

%% === Step 4. 训练逻辑回归模型 ===
glm_model = fitglm(train_data, 'ResponseVar','employment_status', ...
    'PredictorVars', {'age','sex','edu','prof','industry'}, ...
    'Distribution','binomial','Link','logit');
logistic_probs = predict(glm_model, test_data(:,1:5));

%% === Step 5. 训练随机森林模型 ===
rf_model = fitcensemble(X_train, y_train, 'Method','Bag', ...
    'NumLearningCycles',100,'Learners','Tree');
[~, rf_probs] = predict(rf_model, X_test);
rf_probs = rf_probs(:,2);

% 特征重要性图
figure;
bar(oobPermutedPredictorImportance(rf_model));
xticklabels({'age','sex','edu','prof','industry'});
title('随机森林特征重要性');

%% === Step 6. 组合模型预测 ===
% 加权平均（权重可调）
w_log = 0.5; w_rf = 0.5;
combined_probs = w_log * logistic_probs + w_rf * rf_probs;
combined_pred = double(combined_probs > 0.5);

%% === Step 7. 模型评估 ===
function evaluate_model(y_true, y_pred, prob_scores, model_name)
    cm = confusionmat(y_true, y_pred);
    accuracy = sum(diag(cm)) / sum(cm, 'all');
    precision = cm(2,2) / (cm(2,2)+cm(1,2));
    recall = cm(2,2) / (cm(2,2)+cm(2,1));
    f1 = 2 * (precision * recall) / (precision + recall);

    fprintf('\n[%s] 准确率: %.2f%%\n', model_name, accuracy*100);
    fprintf('精确率: %.2f%%, 召回率: %.2f%%, F1: %.2f\n', precision*100, recall*100, f1);

    % 混淆矩阵
    figure;
    confusionchart(cm, {'失业','就业'}, ...
        'Title', sprintf('%s - Confusion Matrix', model_name), ...
        'RowSummary','row-normalized','ColumnSummary','column-normalized');

    % ROC 曲线
    [X,Y,~,AUC] = perfcurve(y_true, prob_scores, 1);
    figure;
    plot(X,Y,'LineWidth',2); grid on;
    xlabel('False Positive Rate'); ylabel('True Positive Rate');
    title(sprintf('%s - ROC (AUC=%.2f)', model_name, AUC));
end

% 分别评估三个模型
evaluate_model(y_test, double(logistic_probs > 0.5), logistic_probs, '逻辑回归');
evaluate_model(y_test, predict(rf_model, X_test), rf_probs, '随机森林');
evaluate_model(y_test, combined_pred, combined_probs, '组合模型');





% 1. 读取预测集数据
predict_data = readtable('预测集.xlsx');

% 2. 年龄段分类
age_bins = [0, 25, 35, 45, Inf];
age_labels = {'<25', '25-35', '35-45', '>45'};
predict_data.age_group = discretize(predict_data.age, age_bins, 'categorical', age_labels);

% 3. 性别转换
sex_labels = {'男', '女'};
predict_data.sex = categorical(predict_data.sex, [1, 2], sex_labels);

% 4. 学历转换
edu_map = containers.Map(...
    {'10','11','12','13','14','15','16','17','18','19', ...
     '20','21','22','23','28', ...
     '30','31','32','33', ...
     '40','41','42','43','44','45','46','47','48','49','91','92','93', ...
     '50','60','61','62','63','70','71','73','80','81','83','90', ...
     '99'}, ...
    {'研究生及以上','研究生及以上','研究生及以上','研究生及以上','研究生及以上',...
     '研究生及以上','研究生及以上','研究生及以上','研究生及以上','研究生及以上',...
     '本科','本科','本科','本科','本科',...
     '专科','专科','专科','专科',...
     '中职','中职','中职','中职','中职','中职','中职','中职','中职','中职','中职','中职','中职',...
     '高中及以下','高中及以下','高中及以下','高中及以下','高中及以下',...
     '高中及以下','高中及以下','高中及以下','高中及以下','高中及以下',...
     '高中及以下','高中及以下',...
     '其他'});
predict_data.edu = cell(height(predict_data), 1);
for i = 1:height(predict_data)
    edu_key = num2str(predict_data.edu_level(i));
    if isKey(edu_map, edu_key)
        predict_data.edu{i} = edu_map(edu_key);
    else
        predict_data.edu{i} = '其他';
    end
end
predict_data.edu = categorical(predict_data.edu);

% 5. 专业简化：统一为训练数据集中的四种类型
predict_data.prof = predict_data.c_aac183;
predict_data.prof = strrep(predict_data.prof, '其他学科', '其他');
prof_categories = unique(predict_data.prof);
if length(prof_categories) > 4
    prof_counts = groupcounts(predict_data, 'prof');
    [~, idx] = sort([prof_counts.GroupCount], 'descend');
    top_profs = prof_categories(idx(1:min(4, length(idx))));
    predict_data.prof(~ismember(predict_data.prof, top_profs)) = {'其他'};
end
predict_data.prof = categorical(predict_data.prof);

% 6. 行业设默认值
predict_data.industry = repmat({'未知'}, height(predict_data), 1);
predict_data.industry = categorical(predict_data.industry);
predict_data.industry(predict_data.industry == '未知') = '服务业及其他';

% 7. 分类变量顺序统一
predict_data.sex = categorical(predict_data.sex, categories(data.sex)); % 与训练集对齐
predict_data.edu = categorical(predict_data.edu, categories(data.edu));
predict_data.prof = categorical(predict_data.prof, categories(data.prof));
predict_data.industry = categorical(predict_data.industry, categories(data.industry));

% 8. 构造特征表用于预测
if ~ismember('people_id', predict_data.Properties.VariableNames)
    % 如果数据缺少 people_id 字段，则添加简单 ID
    predict_data.people_id = arrayfun(@(x) sprintf('ID_%d', x), (1:height(predict_data))', 'UniformOutput', false);
end

% 构造用于预测的特征表
predict_data_table = table(...
    predict_data.age, ...
    double(predict_data.sex), ...
    double(predict_data.edu), ...
    double(predict_data.prof), ...
    double(predict_data.industry), ...
    'VariableNames', {'age','sex','edu','prof','industry'});

% 检查 predict_data_table 中的缺失值
missing_info = ismissing(predict_data_table);  % 获取缺失值逻辑矩阵

% 打印缺失值信息
if any(missing_info, 'all')
    fprintf('预测数据表中存在缺失值，具体信息如下：\n');
    
    % 获取有缺失值的行和列索引
    [rows, cols] = find(missing_info);
    
    % 遍历并打印每个缺失值的具体位置及列名称
    for i = 1:length(rows)
        fprintf('第 %d 行，第 %d 列 (%s) 缺失值\n', ...
                rows(i), cols(i), predict_data_table.Properties.VariableNames{cols(i)});
    end
    
    % 打印含缺失值的完整行，方便调试
    fprintf('\n结合上下文查看带缺失值的完整行：\n');
    disp(predict_data_table(unique(rows), :));
    
    % 填充缺失值：专门针对 'prof' 列 填充为 '其他'
    if ismember('prof', predict_data_table.Properties.VariableNames)
        missing_prof_rows = rows(cols == find(strcmp(predict_data_table.Properties.VariableNames, 'prof')));
        
        % 将 'prof' 列的缺失值填充为 '其他'
        predict_data_table.prof(missing_prof_rows) = categorical({'其他'});
        fprintf('\n缺失值已填充。\n');
    else
        fprintf('\n警告：未找到 "prof" 列，无法填充缺失值。\n');
    end
else
    fprintf('预测数据表中无缺失值。\n');
end

% 9. 使用逻辑回归模型预测概率
prob_logistic = predict(glm_model, predict_data_table);

% 10. 使用随机森林模型预测概率
X_predict = predict_data_table{:,:}; % 转换为数值矩阵
[~, prob_rf] = predict(rf_model, X_predict);

% 11. 模型融合（加权平均）
alpha = 0.5;  % logistic 权重
beta = 0.5;   % random forest 权重
fused_prob = alpha * prob_logistic + beta * prob_rf(:,2);

% 12. 输出预测结果（带人员ID + 分类标签）

threshold = 0.5;

% 12. 输出预测结果（带人员ID + 分类标签）
predict_result = table();
predict_result.people_id = predict_data.people_id;
predict_result.prob_logistic = prob_logistic;
predict_result.prob_rf = prob_rf(:,2);  % 第二列为“就业”概率
predict_result.fused_prob = fused_prob;
predict_result.predicted_label = categorical(fused_prob > threshold, [false true], {'失业', '就业'});

% 显示前几行结果
disp(predict_result(1:10,:));

% 可选：保存为Excel文件
writetable(predict_result, '预测结果输出.xlsx');